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Abstract 


This  document  looks  at  two  techniques  for  modelling  the  radiation  pattern  of  a  parabolic 
reflector  antenna.  Both  these  methods  use  physical  optics  to  determine  the  induced  currents 
on  the  reflector  due  to  the  feed.  The  first  method  uses  MATLAB’s  numerical  integration 
routine  to  compute  the  far-field  radiation  pattern.  The  second  method  involves  segmenting 
the  reflector  surface  and  approximating  the  amplitude  and  phase  of  the  surface  currents  by  a 
first  degree  polynomial.  This  approximation  allows  the  integration  to  be  performed  analyti¬ 
cally  in  closed-form.  Three  different  reflector  configurations  arc  investigated  including  two 
configurations  with  a  laterally  displaced  feed.  Both  these  techniques  yield  accurate  results 
but  the  latter  technique  had  a  significant  improvement  in  computational  efficiency. 


Resume 


Le  present  document  traite  de  deux  methodes  de  modelisation  du  diagramme  de  rayonne- 
ment  d’ une  antenne  a  reflecteur  parabolique.  Les  deux  methodes  font  appel  a  l’optique  phy¬ 
sique  pour  determiner  les  courants  induits  sur  le  reflecteur  par  T  alimentation.  La  premiere 
methode  utilise  une  routine  d’integration  numerique  MATLAB  pour  calculer  le  diagramme 
de  rayonnement  en  champ  lointain.  La  deuxieme  effectue  une  segmentation  de  la  surface 
du  reflecteur  et  une  approximation  de  1' amplitude  et  de  la  phase  des  courants  de  surface 
au  moyen  d’un  polynome  du  premier  degre.  Cette  approximation  permet  d’effectuer  une 
integration  analytique  en  forme  fermee.  Trois  configurations  de  reflecteur  differentes  sont 
etudiees,  y  compris  deux  configurations  avec  une  alimentation  deplacee  lateralement.  Ces 
deux  methodes  donnent  des  resultats  precis,  mais  la  deuxieme  presente  une  efficacite  com- 
putationnelle  nettement  superieure. 
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Executive  summary 


DRDC  Ottawa  is  investigating  the  design  of  missile  models  that  can  be  used  to  predict 
missile  performance  in  the  presence  and  absence  of  Electronic  Countermeasures  (ECM). 
This  is  the  first  of  a  series  of  publications  that  will  deal  with  the  different  components  that 
make  up  the  seeker  and  guidance  system  of  a  radar  guided  anti-ship  missile.  This  publi¬ 
cation  examines  different  methods  of  modelling  the  missile  seeker  antenna  accurately  and 
efficiently.  Once  an  appropriate  antenna  modelling  method  has  been  identified,  it  will  be 
integrated  into  a  larger  simulation  framework  that  attempts  to  capture  the  operation  of  the 
complete  missile  in  its  operating  environment.  While  antenna  modelling  is  computation¬ 
ally  intense  it  is  important  that  the  chosen  model  not  significantly  slow  down  the  missile 
engagement  simulation.  Determining  the  best  ECM  technique  requires  a  large  number  of 
simulation  runs,  therefore  it  is  important  that  these  individual  runs  execute  in  seconds  ver¬ 
sus  hours. 

This  document  looks  at  two  techniques  for  modelling  the  radiation  pattern  of  a  parabolic 
reflector  antenna.  Both  these  methods  use  physical  optics  to  determine  the  induced  currents 
on  the  reflector  due  to  the  feed.  The  first  method  uses  MATLAB’s  numerical  integration 
routine  to  compute  the  far- Held  radiation  pattern.  The  second  method  involves  segmenting 
the  reflector  surface  and  approximating  the  amplitude  and  phase  of  the  surface  currents  by 
a  first  degree  polynomial.  This  approximation  allows  the  integration  to  be  performed  ana¬ 
lytically  in  closed-  form.  Three  different  reflector  configurations  arc  investigated  including 
two  configurations  with  a  laterally  displaced  feed.  Both  of  the  methods  discussed  (Physical 
Optics/Surface  Current  (PO/SC)  Method  and  Ludwig’s  Method)  give  accurate  results  when 
compared  to  patterns  computed  using  commercially  available  reflector  modelling  software. 
Where  these  two  methods  differ  is  in  the  amount  of  time  required  to  compute  these  re¬ 
sults.  The  PO/SC  Method,  with  a  lateral  displacement  of  the  feed,  where  the  integrand  of 
the  radiation  integral  increases  in  complexity,  requires  in  excess  of  one  hour  to  evaluate  one 
principal  plane  pattern.  Incorporating  this  model  into  a  seeker  model  would  therefore  be  im¬ 
practical.  As  an  alternative,  Ludwig’s  Method,  regardless  of  the  lateral  displacement  of  the 
feed,  can  accurately  compute  the  principal  plane  pattern  in  less  than  10  seconds.  Because  of 
its  speed  and  accuracy,  future  work  will  involve  applying  Ludwig’s  method  to  monopulse 
reflector  antennas  having  a  compound  feed.  Finally  this  more  efficient  monopulse  antenna 
model  could  then  be  incorporated  into  a  missile  seeker  model. 

A  number  of  reports  will  be  written  over  the  next  year  that  will  provide  the  foundation 
for  future  countermeasures  research.  These  reports  are  the  result  of  studies  carried  out  as 
part  of  Project  llat  and  will  be  used  to  develop  future  Thrust  11a  Projects  that  will  be 
required  to  provide  support  for  the  development  of  EW  techniques  in  the  future.  By  pro¬ 
viding  a  deeper  understanding  of  the  interactions  between  the  jammer  and  threats,  these 
studies  will  also  help  to  define  the  requirements  for  future  Naval  EW  systems. 
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Sommaire 


RDDC  Ottawa  mene  des  recherches  sur  la  conception  de  modeles  de  missile  susceptibles 
de  servir  a  la  prediction  du  rendement  des  missiles  en  presence  et  en  1’ absence  de  contre- 
mesures  electroniques  (CME).  Le  present  document  est  le  premier  d’une  serie  de  publi¬ 
cations  sur  les  differents  composants  de  l’autodirecteur  et  du  systeme  de  guidage  d’un 
missile  antinavires  guide  au  radar.  La  presente  publication  examine  differentes  methodes 
de  modelisation  de  l’antenne  de  l’autodirecteur  du  missile  avec  precision  et  de  facon  effi- 
ciente.  Une  fois  qu’une  methode  appropriee  de  modelisation  de  l’antenne  aura  ete  identifiee, 
elle  sera  integree  a  un  cadre  de  simulation  plus  vaste  qui  saisit  autant  que  possible  le  fonc- 
tionnement  du  missile  complet  dans  son  milieu  operationnel.  Bien  que  la  modelisation  de 
l’antenne  soit  intense  sur  le  plan  des  calculs,  il  est  important  que  le  modele  choisi  ne  ra- 
lentisse  pas  significativement  la  simulation  de  l’engagement  du  missile.  Pour  determiner  la 
meilleure  technique  CME,  il  faut  un  grand  nombre  de  passages  en  simulation ;  c’est  pour- 
quoi  il  est  important  que  les  passages  individuels  soient  executes  en  quelques  secondes 
plutot  qu’en  quelques  heures. 

Le  present  document  traite  de  deux  methodes  de  modelisation  du  diagranmie  de  rayonne- 
ment  d’ une  antenne  a  reflecteur  paraboliquc.  Les  deux  methodes  font  appel  a  l’optique  phy¬ 
sique  pour  determiner  les  courants  induits  sur  le  reflecteur  par  L  alimentation.  La  premiere 
methode  utilise  une  routine  d’integration  numerique  MATLAB  pour  calculer  le  diagranmie 
de  rayonnement  en  champ  lointain.  La  deuxieme  effectue  une  segmentation  de  la  surface 
du  reflecteur  et  une  approximation  de  1' amplitude  et  de  la  phase  des  courants  de  surface 
au  moyen  d’un  polynome  de  premier  degre.  Cette  approximation  permet  d’effectuer  une 
integration  analytique  en  forme  fermee.  Trois  configurations  de  reflecteur  differentes  sont 
etudiees,  y  compris  deux  configurations  avec  une  alimentation  deplacee  lateralement.  Les 
deux  methodes  examinees  (methode  de  determination  du  courant  de  surface  a  l’aide  de 
l’optique  physique  et  methode  de  Ludwig)  donnent  des  resultats  precis  par  rapport  aux  dia- 
grammes  calcules  au  moyen  de  logiciels  commerciaux  de  modelisation  du  reflecteur.  La 
ou  ces  deux  methodes  different,  c’est  dans  le  temps  requis  pour  le  calcul  des  resultats.  La 
methode  de  determination  du  courant  de  surface  a  l’aide  de  l’optique  physique,  qui  tient 
compte  du  deplacement  lateral  de  1’ alimentation,  ou  la  fonction  a  integrer  de  l'integrale 
de  rayonnement  augmente  en  complexity  requiert  plus  d’une  heure  pour  l’e valuation  d’un 
diagramme  sur  le  plan  principal.  Il  ne  serait  done  pas  pratique  d'integrer  cette  methode  a 
un  modele  d’autodirecteur.  Par  contre,  la  methode  de  Ludwig,  peu  importe  le  deplacement 
lateral  de  1’ alimentation,  peut  calculer  avec  precision  le  diagramme  sur  le  plan  principal  en 
moins  de  10  secondes.  Grace  a  la  rapidite  et  a  la  precision  de  cette  methode,  les  recherches 
futures  feront  appel  a  1’ application  de  techniques  de  modelisation,  fondees  sur  la  methode 
de  Ludwig,  aux  antennes  a  reflecteur  ayant  une  alimentation  composee  pour  utilisation  dans 
des  systemes  monopulses.  Enfin,  un  modele  d’ antenne  monopulse  pourrait  alors  etre  integre 
a  un  modele  d’autodirecteur  de  missile. 

Au  cours  de  l’annee  qui  vient,  un  certain  nombre  de  rapports  seront  rediges  pour  mettte  en 
place  les  fondements  des  recherches  futures  sur  les  contre-mesures.  Ces  rapports  seront  les 
resultats  d’etudes  menees  dans  le  cadre  du  projet  1  lat  et  serviront  a  la  definition  des  projets 
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futurs  du  volet  1  la  qui  seront  requis  a  l’appui  de  la  mise  au  point  de  techniques  de  guerre 
electronique  (GE)  a  l’avenir.  En  permettant  une  meilleure  comprehension  des  interactions 
entre  le  brouilleur  et  les  menaces,  ces  etudes  contribueront  aussi  a  la  definition  des  besoins 
en  ce  qui  concerne  les  futurs  systemes  GE  navals. 


Frederic  Arpin.  2005.  Parabolic  Reflector  Modelling  Techniques  with  a  Laterally  Displaced 
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1  Introduction 


This  document  will  investigate  two  methods  for  the  modelling  of  the  radiation  pattern  of  a 
parabolic  reflector  with  a  laterally  displaced  feed.  Both  techniques  arc  based  on  physical 
optics  and  both  techniques  involve  determining  the  surface  currents  induced  on  the  reflector 
surface. 

The  first  technique  will  integrate  the  induced  surface  currents  over  the  reflector  surface 
which  will  yield  the  far-field  radiation  pattern.  This  integration  will  be  performed  us¬ 
ing  MATLAB's  built-in  numerical  integration  routine  that  uses  recursive  adaptive  Simpson 
quadrature  to  approximate  the  integral. 

The  second  method,  instead  of  integrating  the  surface  currents  directly,  imposes  a  inte¬ 
gration  grid  on  the  surface  of  the  reflector,  where  within  each  grid  element  the  amplitude 
and  phase  of  the  integrand  is  approximated  by  a  simple  linear  function.  This  function  over 
the  incremental  surface  area,  which  make  up  a  grid  element,  can  be  integrated  analytically. 

The  goal  here  is  to  determine  a  computationally  efficient  modelling  technique  that  can  be 
used  to  model  monopulse  tracking  antennas  which  can  then  be  incorporated  into  a  missile 
seeker  model. 

2  Modelling  of  Far-Field  Radiation  Patterns  for 
Reflector  Antennas 
2.1  PO/Surface  Current  Method  (PO/SC) 

Consider  the  problem  geometry  shown  in  Figure  1.  To  determine  the  radiation  pattern  of 
this  reflector  we  make  use  of  the  far-field  radiation  integral  given  as  [1] 

Erad(epAp)  =  -jkzf  '^  (I- rprp)  J  Jj^rds  (1) 

Where  fp,  rs  arc  defined  in  Figure  1,  k  is  the  free  space  propagation  constant,  Z0  is  the 
intrinsic  impedance  of  free  space,  7  is  a  unit  dyadic  and  Js  is  the  induced  surface  current  on 
S. 

In  order  to  evaluate  (1)  we  must  first  determine  the  surface  current  Js,  induced  onto  the 
reflector  surface,  by  the  means  of  the  physical  optics  (PO)  approximation.  The  PO  approx¬ 
imation  makes  three  assumptions,  namely: 

•  The  reflector  has  an  electrically  large  radius  of  curvature  so  that  locally  at  each  reflec¬ 
tion  point  the  surface  of  the  reflector  can  be  considered  planar. 

•  The  radius  of  curvature  of  the  incident  field  is  large  and  locally  at  each  reflection  point 
can  be  viewed  as  a  plane  wave. 
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Asi'r? 


Figure  1:  Parabolic  Reflector  Problem  Geometry 


•  The  reflector  is  considered  to  be  a  perfectly  conducting  surface  which  yields  a  reflected 
field  equal  to  the  incident  field. 

If  we  then  have  a  locally  plane  wave  incident  on  a  locally  plane  perfectly  conducting  sur¬ 
face,  referring  to  Figure  2,  the  boundary  conditions  tell  us  that  the  current  induced  on  the 


(7  =  °° 

Figure  2:  Induced  Currents  on  a  Perfect  Conductor 


surface  of  the  conductor  will  be  given  by 

J~s  =  Htan  =fixH  =  nx  (H  +  Hr)  (2) 

Now  considering  the  third  assumption  above  we  have  that  the  reflected  field  is  equal  to  the 
incident  field,  we  can  then  write  the  induced  surface  current  as 

Js  =  n  x  (Hi  +  Hj)  =  2 h  x  H  (3) 

We  can  now  proceed  to  determining  the  equations  defining  Js.  We  will  start  by  defining  the 
electric  field  from  the  feed  antenna,  for  the  purpose  of  this  report  we  will  model  the  feed  in 
terms  of  its  principal  plane  patterns.  For  a  y-polarized  feed  we  have 

_ .  _  „  „ 

Emc(?  s)  =  [(c£(e;,)Si<)e'+(ctf(e:)Co4;)fs]  —  (4) 

rs 
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where  Ce(0')  and  0/(0')  define  the  E-Plane  and  H-Plane  principal  patterns  respectively, 
which  are  defined  as 

C£(0'v)  =  {cosQ'syE 

CM)  =  {cos$yH 

qE  and  qH  are  selected  such  as  to  match  the  CE  ( 0( )  and  Cf/fO'J  to  the  principal  plane 
patterns  of  the  actual  feed.  The  magnetic  incident  field  can  be  related  to  the  incident  electric 
field  by  the  following  equation 


Hinc{Ps) 


^Ps  x  Einc(Ps) 


(5) 


We  can  then  determine  the  expression  for  the  incident  magnetic  field  given  (4)  as  follows 


?s  x  Einc(Ps)  = 


p  s  0',  <K, 

1  0  0 
0  Ee;  Eys 


1 


,-jkp 


Hinc{?s)  =  -PsxEi"\Ps)=[{-CH{MosMs  +  {CE(MMsW.^  7  , 

^ O V s 

This  can  now  be  expressed  in  cartesian  components  as 


(6) 


Hmc(Ps )  =  \(—CHCos%'scos2§'s  —  CEsin2§'s)xs 

+  ( — CHCOsQ'scosty'ssin§'s  +  CECOs§'ssin§'s)y's 

e~Jkr * 

+  (CHsinQ'scos^s)zs] 

^ o's 


(V) 


We  will  next  define  the  unit  normal  h  to  the  reflector  surface.  We  are  considering  a  parabolic 
reflector  and  the  surface  is  defined  as 


2  F 

1  +  CO.V0* 


(8) 


where  F  is  the  focal  length  of  the  parabola.  (8)  can  be  rewritten  as 


0  =  2  F  —  rs  —  rscosQs 


(9) 


The  surface  of  the  reflector  is  then  defined  as 


5  =  2/  —  r,  —  rscosQs 

The  normal  of  this  surface  can  be  determined  by  taking  the  gradient  of  5  as  follows 

55  „  1  55  A  1  55, 

n  —  VS  —  —rs  -| —  ^-0S  H - r^r 


we  then  have 


5 rs  ‘  rs  50,  rssinQs  5(), 


n  =  —  (1  +  cosQs)rs  +  (sinds)ds 


(10) 


(11) 

(12) 
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Now  the  unit  normal  is  given  by 


n  = 


We  then  have,  using  several  trigonometric  identities  and  simplifications 

n  =  (—sin(Qs/2)costys)xs  —  (, sin(Bs/2)sintys)ys  —  (cos(Qs/2))zs 


(13) 


(14) 


which  is  the  normal  that  is  required  to  determine  the  induced  current  given  in  (3).  It  should 
be  noted  that  xs  =  x's,  ys  =  y's  and  zs  =  z's  which  is  evident  in  Figure  1.  Using  (eq2. 1.7)  and 
(eq2.1.3)  in  (eq2.1.3)  we  have 


h  x  Hinc{? s)  = 


xs  }’s  zs 


Hx  Hy  Hz 


n  x  Hinc{?s)  =  ( nyHz  -  nzHy)xs  +  (nzHx  -  nxHz)ys  +  (nxHy  -  nyHx)zs  (15) 

(nyHz  —  nzHy)  =  [ — C//.sm  ( 9S  /  2.)sin§ssinQ'scos§'s 
—  Chcos{Qs /2)cos$scos§'ssin§'s 

e~Jk>Js 


+  Cecos  (0,/2)  cos§'ssin§'s 


Z0r's 


(16) 


(nzHx  —  nxHz)  =  \Chcos{Qs/2)cos2§'scos$s 
+  C[7  cos  ( 0 ,  /2 )  sin2$'s 

+  Cj]sin{Qs/2.)cos§scos§'ssin§'^\ 


e-jkK 

7  r' 

^O'  s 


(17) 


(nxHy  —  nyHx)  =  \CHsin(Qs/2)cos§scosQlscos§'ssin§'s 
—  Cesui{Qs  /2)cos§scos§'ssin§'s 
—  CHsin(Qs/2)sin§scos2§'scosQ's 

e-jkK 


2  A/1 


—  C[;sin(Qs/2)sintyssinzty 


7  r' 

^o'  s 


We  then  have 


f(x,y,z)  =  2h  x  Hmc(?s) 


Therefore, 


=  2 (nyHz  -  nzHy)xs  +  2 (nzHx  -  nxH-)ys  +  2 (nxHy  -  nyHx)zs 


J x  =  [— CHsin(Qs/2)sin§ssinQ'scos§'s 
—  Chcos{Qs /2)cos$scos§'ssin§'s 

2e~ikr>s 

+  CEcos(Qs/2)costy'ssirn |>']  - 

70rs 


(18) 


(19) 


(20) 
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Jy  =  \Chcos(Qs/2)cos2§'scosQ's 
+  CEcos{Qs/2)sin2^s 

+  CHsin(Qs  /  2)cos§scos§,ssin§'s\ 


(21) 


Jz  =  [CHSin(Qs/2)cOS§sCOsQ'sCOS§'ssin§'s 
—  CEsin  (0,  /  2)  co^cost))'  sin§'s 
—  Cfisin(Qs/2)  sin§scos 2  c^cas'd^. 

2  e~ikf 

-  CE sin (0,v/2) sin§ssin2§'s]  -  (22) 

Z„ts 

Now  prior  to  integrating  these  induced  surface  currents  we  must  define  the  primed  quanti¬ 
ties,  namely  r's,  0'  and  (|)',  in  terms  of  rs,  0S,  (f>v  which  are  the  known  quantities  based  on  the 
reflector  geometry.  The  vector  r'iV  is  given  by 

r's  =  rs-f  (23) 

Note  that  /  is  the  vector  that  defines  the  feed  position  and  is  given  as 

f  =  xfxs+yfys  (24) 


and  rs  is  given  as 


rs  =  ( rssinQscostys)xs  +  (rssinQssin$s)ys  +  ( rscosds)zs 


(25) 


We  then  have 


r's  =  (rssinQscosfys  -  fx)xs  +  (rssin8ssin§s  -  fy)ys  +  (. rscosQs)zs 

Using  (26)  we  can  now  define  /s,  0'  and  (])'  as  follows 


and 


r's=  \/r'2sx  +  r'2y  +  tj2z 


0;v  =  tan 


'rs+fx+fy  ~  2rssinQs(fxcos<\)s  +  fysin§s) 

0;v  =  tan  1 

\j  r'lx  + 

y! 

'  sz 

|"  sj  (rssinQscos§< 

:  -  fx)2  +  (rssinQssintys  -  fy)2 

t\cosQ v 


^=tan  1  (i^) 


c[)j  =  tan 


!  f  rssinQssinfys  —  fy 
rssinQscos§s  -  fx 


(26) 


(27) 


(28) 


(29) 
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Prior  to  integrating  the  induced  surface  currents  (20)-(22)  we  also  need  to  define  the  two 
additional  terms  in  the  integration  shown  in  equation  (1),  namely  the  exponential  involving 
a  dot  product  and  the  differential  surface  ds. 

First  let  us  define  the  dot  product  in  the  exponential  term,  we  have 

F  '  p  —  >'sx ^ px  T  FsyFpy  T  ^sz^pz  (30) 


where 


rs  =  (. rssinQscostys)xs  +  (rssinQssinfys)ys  +  (rscosQs)zs 


and 

fp  =  ( sinQpcosfyp)xp  +  ( sinQpsintyp)yp  +  (cosQp)zp 

Before  performing  the  dot  product  we  must  first  express  rp  in  the  xs,ys, zs  coordinate  system, 
which  gives 

fp  =  ( sinQpcos§p)xs  —  (sinQpsin§p)ys  —  (cosQp)zs 


We  then  have 


rs  ■  fp  =  rssinQscostyssinQpcostyp  —  rssinQssinfyssinQpsintyp  —  rscosQscosQp  (31) 
Next  the  differential  surface  area  ds  is  defined  as 

ds  =  r^sec(Qs/2)sin(Qs)dQsd§s  (32) 


Expressions  (20-22)  along  with  (31)  and  (32)  can  now  be  numerically  integrated  separately 
yielding 


Fx  =  1  jj^'+'ds 

(33) 

Fy=  1  jjyejkT^’ds 

(34) 

FZ  =  J  j^e^nls 

(35) 

Now  that  the  integral  has  been  evaluated  we  can  finish  evaluating  expression(l).  To  do  this 
we  first  must  evaluate  the  following: 

d-rprp)-F 

(36) 

which  can  be  rewritten  as 

F-rP(rP-F) 

(37) 

It  should  be  noted  that  F  is  defined  in  the  xs,ys,zs  coordinate  system  and  rp 
the  xp,yp,zp  system,  as  a  consequence  rp  in  the  xs,ys,zs  system  is  given  as 

is  defined  in 

fp  =  (sinQpcosfyp)xs  —  (sinQpsinfyp)ys  —  (cosdp)zs 

(38) 

We  then  have 

fp  F  =  (Fxsin9pcos§p)  +  (— FysinQpsin§p )  +  (—  FycosQp) 

(39) 
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Then 


fp(fp  ■  F)  =  ( Fxsin2Qpcos2§p  —  Fysin2Qpsintypcostyp  —  FzcosQpsinQpcos(. \>p)xp 
+  (Fxsin2QpcostypSinfyp  —  Fysin2Qpsin2typ  —  FzcosQ  psinQ  psin<\>  p)y  p 
+  ( FxsinQpcosQpcostyp  —  FysinQpcosQpsintyp  —  Fzcos2Qp)z.p  (40) 

We  can  now  evaluate  F  —  fp(fp-F)  however  F  is  defined  in  the  xs,ys,  zs  system  and  rp (rp  ■  F ) 
is  defined  in  xp,yp,zp  system.  We  can  then  define  F  in  the  xp,yp,zp  system  and  we  have 
F  =  Fxxp  —  Fyyp  —  Fzzp,  we  then  have 

F  —  fp(fp  •  F)  =  (Fx  —  Fxsin2QpCos2§p  —  Fysin2Qpsin§pcos§p  —  FzcosQpsinQpcos§p)xp 

—  (Fy  +  Fxsin2Qpcos<\)pSin<\)p  —  Fysin2Qpsin2<\ \>p  —  FzcosQ  psinQ  psin<\>  p)y  p 

—  (Fz  +  FxsinQ  pcosQ  pcosfy  p  —  Fy  s  in  0  p  co  s  0  p  s  i  n§  p  —Fzcos2  Qp)zP  (41) 

Finally  the  reflector  pattern  is  given  as 


II 

1 

Jr- 

3  4.. 

13* 

[(F-rp(fp-F))-xp] 

e~jkrp 

II 

1 

c 

3 

[(F-rp{rp-F))-yp\ 

II 

1 

^  1 

tV5  *-"•  1 

IS8 

[(F-rp(rp-F))-zP\ 

(42) 

The  radiation  pattern  can  also  be  expressed  in  its  co-pol  and  x-pol  components  as  follows  [2] 

E™d(Qp^p)  =  [-(1  -cosQp)sin§pcos§p]Exad 
+  [1  —sin2§p{\  —cosQp)]  Eyad 
—  [sinQpsin§p\Ezad  (43) 

Exad(Qp,§p)  =  [1  -  cos\(  1  -  cosQp)]Exad 

—  [(1  —  cosQp)sintypcosfyp]E'yad 

—  [sinQpCosfyp]  Ezad  (44) 

2.2  Ludwig’s  Method 

The  Ludwig  method  is  one  of  the  first  numerical  techniques  to  evaluate  the  radiation  inte¬ 
gral,  with  the  first  work  dating  back  to  1968  [3].  Ludwig’s  method  also  makes  use  of  the 
PO  approximation  and  its  stalling  point  is  the  radiation  integral  given  in  equation  (1)  as 
follows 

Erad(Qp^p)  =  -jkZ0— -(f-rprp)  [  [ Jsejkf^ds  (45) 

r p  J  Js 

The  integration  of  equation  (45)  can  be  rewritten  in  a  more  general  form  as 

F(QPAP)  =  j  (46) 
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S"(0s  ,(])y)  is  given  by  equations  (20)-(22),  from  the  previous  section,  with  the  exponential 
terms  removed  and  the  non-differential  portion  of  equation  (32).  The  exponential  term  will 
be  lumped  with  the  y(0v,(|)v,0p,(|)/,)  term  of  equation  (46).  We  then  have  K(QS,§S),  in  its 
three  cartesian  components  given  as 


kx{  e.v,^v)  = 


+ 


[— Chsui  (0lS /2)  sin§ssinQ'scos§'s 
CHCOs(Qs/2)cosQ,scos§'ssin§'s 

C£«M(e!/2)C»<Sm<J^“C(e'/2)s,'"(e’ 


7  r' 
^o'  s 


(47) 


Ky{9s,§s)  = 


+ 

+ 


\Chcos(Qs/2)cos2§'scosQ's 

CECOs($s/2)sin2§'s 


CHsin(Qs/2)cos§scos§/ssinQ[. 


2r2sec(Qs /2)sin(Qs) 


(48) 


A7(0iV,(|)y)  =  [CHsin(Qs/2)cos§scosQ'scos§'ssin§'s 
—  CEsin(Qs/2)cos§scos§/ssin§/s 
—  Cnsin(Qs /2)sin§scos2§'s,cosQ[, 


—  CEsin(Qs/2)sin<\>ssin<\)s 


2,n  2r2sec(Qs/2)sin(Qs 


7  r' 

^O'  S 


(49) 


and  y(Qs,§s,Qp,§p)  will  be  given  by 


Y(0i > 4*^ 5  —  rs  4"  rs '  t'p 


(50) 


rs  ■  fp  is  given  by  equation  (31)  of  the  previous  section. 


Now  consider  a  fixed  output  radiation  pattern  point  (Qpp,§qp)  for  a  single  cartesian  com¬ 
ponent  of  F  and  K,  the  integration  is  then 

/•e?  rtf 

Fpq=  /  K(QMeJ^^dQd^  (51) 

Jo  Jo 

Here  K  is  specified  on  an  M  x  N  integration  grid  for  the  points  (0'”,(|)").  We  can  now 
consider  how  the  integrand,  (K(QS,<. behaves  over  an  incremental  surface  A  Smn 
on  the  surface  of  the  reflector  where 


A Smn  =  { (©„<!>,)  :  0'f  <  0.V  <  0"'+1 ,  C  <  ^  <  C+1 }  (52) 

If  A  Smn  is  physically  on  the  order  of  a  wavelength  the  electrical  path  length  term  (y/cy(0s.(|)v)) 
cannot  vary  by  more  than  2k.  Now  consider  the  behavior  of  K(QS.<\>S)  over  the  incremental 
area  A Smn.  Since  typically  electromagnetic  fields  do  not  have  abrupt  changes  over  a  dis¬ 
tance  on  the  order  of  a  wavelength,  we  can  then  say  that  the  integrand,  (K(QS,(. tys)e-'kv(®s^s>). 
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will  not  vary  abruptly  and  be  well  behaved  over  the  incremental  area  A Smn.  This  being  said 
K  and  yean  be  separately  approximated  by  a  linear  function  over  A  Smn.  We  then  have 

K(QSi  C)  —  & mn  T  bmn  ( 9V  9  V  )  T  Cmn  (tjh  C ) 

Y(0*,  O  ^  a mn  +  Pmn (Q.v  -  6”')  +  \ mn (C  -  C) 

for(Qs,tys)  e  A Smn  (53) 

The  method  used  in  [3],  for  determining  the  coefficients  of  (53)  is  to  use  a  best  fit  mean- 
squared  plane  to  the  values  of  the  functions,  (K( Qs,§s)e-ik^ds^s'>),  at  the  corners  of  A Smn. 
This  gives 

amn  =  ^  [3K(0™,$)-K(G”,+1,$+1) 

+  ^(e'"+1,C)+^(0”'5C+1)] 

bmn  = 

+  ^(e",+1,C+1)-^(0”'>C+1)] 

Cmn  =  2^[W,C+1)-m",C) 

+  /r(0"'+1,c+1)-^(er+1,  €)]  (54) 

and 

ttrnn  =  ^[3y(0'",C)-Y(0r+1,€+1) 

P  mn  =  ^m+\W-V(e?,W 

+  y(0”,+ 1 1 C+ 1 )  —  y(0”'  >  C+ 1 )  ] 

U  =  ^  [y(0M+Vy(CC) 

+  Y(er+1,C+1)-Y(er+1,C)]  (55) 

where 

A6'"  =  e';i+1-e'" 

AC  =  C+1-C  (56) 


With  (53)-(56)  we  can  now  perform  the  integration  in  closed  form,  and  over  the  incremental 
surface  area  A Smn  we  get  the  following  contribution 


A  F  — 

LA1  mn  — 


jkftmn 


+  b„ 


+  Cn 


A0' 


s_  jk$mnAQ™  _ 


jkfymn 

ejk$m„AQ"'  _  l 
jTP/H« 


jk^>mn 

ejk$mn  A8”  _  y 

(,/T  fimn )  ~ 


_AC 


jk^mn 


(57) 


DRDC  Ottawa  TM  2005-098 


9 


To  get  the  value  of  the  integral  we  simply  need  to  sum  all  these  contributions  over  m  and  n. 
This  can  then  be  repeated  for  each  radiation  pattern  observation  point.  Therefore  in  order 
to  obtain  the  complete  desired  radiation  pattern  we  would  then  have  to  perform  M  xN  x 
P  x  Q  computations.  However  it  should  be  noted  that  equations  (47)-(49),  defining  K  are 
independent  of  the  output  observation  angles  Qp  and  (f>p,  therefore  amn,  bmn  and  cmn  only 
need  to  be  computed  once.  This  will  significantly  reduce  the  computation  time,  however 
it  does  require  some  memory  allocation  for  storage  of  these  coefficients.  We  can  now  use 
equations  (42)-(44),  from  the  previous  section  to  compute  the  complete  co-pol  and  x-pol 
radiation  patterns. 


3  Results 


Each  of  the  formulations  of  the  previous  section  were  coded  in  MATLAB.  We  will  now 
look  at  radiation  patterns  for  parabolic  reflector  antennas  using  these  modelling  methods. 
We  will  specifically  examine  the  co-pol  and  x-pol  patterns  in  the  two  principal  planes  for 
three  different  configurations,  namely 

•  Reflector  #1 

-  D  =  201 

-  F/D  =  0.8 

-  Feed  Location:  (x,y,z)  =  (0,0,0) 

•  Reflector  #2 

-  D  =  201 

-  F/D  =  0.8 

-  Feed  Location:  (x,y,z)  =  (X/4, 0,0) 

•  Reflector  #3 

-  D  =  201 

-  F/D  =  0.8 

-  Feed  Location:  (x,y,z)  =  (0,X/2,0) 

The  radiation  patterns  for  each  of  the  reflector  configurations,  and  each  of  the  methods  de¬ 
scribed  in  the  previous  section,  will  be  compared  to  the  results  from  a  commercial  reflector 
modelling  tool,  GRASP  by  Ticra  [4],  We  will  be  looking  at  the  radiation  patterns  for  301 
points  between  —30°  <  Qp  <  30° 

3.1  Feed  Pattern 

For  the  puipose  of  this  modelling  the  feed  pattern  is  approximated  by  equation  (4),  recall 

Emc(ps )  =  [(c£(0',)™4:)0'+(cff(e:)H:)fj  m 

rs 
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This  approximate  feed  pattern  will  make  use  of  the  principal  plane  patterns  (E-plane  — > 
<]>'  =  0°  and  H-plane  — ►  (j)'  =  90°)  where  Ce(0')  =  cos(Q's)qE  and  0/(0')  =  cos(Q's)qH  rep¬ 
resent  the  E  and  H-plane  patterns  respectively.  We  must  then  select  qE  and  qH  in  order 
to  match  the  principal  plane  patterns  of  the  actual  feed  that  we  are  modelling.  This  feed 
approximation  will  model  the  field  from  the  feed,  for  any  angle  (|)'  by  means  of  interpola¬ 
tion.  As  a  general  rule  [5],  to  achieve  peak  aperture  efficiency,  the  edge  illumination  should 
be  approximately  -lldB.  We  will  then  select  qE  and  qH  that  would  give  us  this  desired 
edge  illumination.  For  the  puipose  of  this  modelling  we  will  look  at  a  balanced  feed  where 
qE  =  qH.  This  type  of  balanced  feed  has  a  very  low  cross-polarized  field  [6]. 

For  the  reflector  we  are  considering  (i.e.  Reflector  #1),  we  will  have  qE  =  qH  =  6.5  and 
with  this  we  have  the  feed  pattern  shown  in  Figure  3.  Note  that  0„  is  the  reflector  sub¬ 


tended  angle,  depicted  in  Figure  4.  The  subtended  angle  is  calculated  using  the  following 
expression 

0O  =  2  tan~ 1  ;  (59) 

4  F 

In  the  case  of  Reflector  #1  the  subtended  angle  0O  =  34.708.  We  can  then  see  from  Figure  3 
that  the  edge  illumination  is  approximately  —11  dB.  It  should  also  be  noted  that  for  the  other 
two  reflector  configuration,  with  the  laterally  displaced  feeds,  we  will  not  be  modifying  the 
feed  model  (i.e.  qE  and  qH  will  remain  unchanged).  As  a  result  we  will  no  longer  have  a 
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Figure  4:  Reflector  Geometry  with  Subtended  Angle  9, 


uniform  —11  dB  edge  illumination  on  the  reflector,  because  of  the  feed  displacement  from 
the  center  axis. 
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3.2  PO/Surface  Current  Method 


This  method,  as  stated  previously,  is  based  on  numerically  integrating  the  surface  currents 
induced  on  the  reflector.  In  order  to  perform  these  integrations  we  make  use  of  MAT- 
LAB ’s  internal  integration  routines,  namely  dlbquad  which  uses  a  recursive  adaptive  Simp¬ 
son  quadrature  to  numerically  evaluate  the  desired  integral.  This  MATLAB  routine  will 
approximate  the  integral  within  a  defined  absolute  error  tolerance  which  by  default  is  set  to 
1  x  10  6.  This  error  tolerance  governs  the  speed  and  accuracy  of  the  numerical  integration. 
A  large  error  tolerance  will  result  in  rapid  computation  with  decreased  accuracy  and  a  small 
error  tolerance  will  provide  very  accurate  results  at  the  expense  of  evaluation  speed.  In  this 
section  we  will  look  how  this  error  tolerance  affects  the  accuracy  and  computing  speed  of 
the  radiation  patterns. 

3.2.1  Reflector  #1 

We  will  first  look  at  how  the  error  tolerance  will  affect  the  accuracy  of  the  radiation  pattern. 
In  Figures  5-9  we  see  the  co-pol  pattern  for  a  single  pattern  cut  (§p  =  90°)  for  5  different  er¬ 
ror  tolerances.  The  computed  patterns  are  compared  to  the  pattern  computed  using  GRASP. 
We  can  clearly  see  from  Figures  5-9  how  the  error  tolerance  affects  the  accuracy  of  the 


Figure  5:  Radiation  Pattern  for  Reflector  #1 ,  MATLAB  Numerical  Integration,  Error 
Tolerance  =  1  x  1(T4,  §p  =  90° 


results.  An  error  tolerance  of  1  x  10  8  was  requires  to  accurately  model  the  radiation  pat¬ 
tern.  Table  1  shows  the  CPU  time  (3.0GFlz  Pentium  IV,  512MB  RAM)  required  to  compute 
each  of  the  patterns  for  the  different  error  tolerances.  We  can  also  see  in  this  table  the  value 
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-70 


-80 - 1 - 1 - 1 - 1 - 1 - 

-30  -20  -10  0  10  20  30 

Angle  off  Broadside  -  0  (degrees) 

Figure  6:  Radiation  Pattern  for  Reflector  #1 ,  MATLAB  Numerical  Integration,  Error 
Tolerance  =  1  x  10-5,  4>p  =  90° 


Figure  7:  Radiation  Pattern  for  Reflector  #1 ,  MATLAB  Numerical  Integration,  Error 
Tolerance  =  1  x  10~6,  §p  =  90° 
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-70 


-80 - 1 - 1 - 1 - 1 - 1 - 

-30  -20  -10  0  10  20  30 

Angle  off  Broadside  -  0  (degrees) 

Figure  8:  Radiation  Pattern  for  Reflector  #1 ,  MATLAB  Numerical  Integration,  Error 
Tolerance  =  1  x  10-7,  4>p  =  90° 


Figure  9:  Radiation  Pattern  for  Reflector  #1 ,  MATLAB  Numerical  Integration,  Error 
Tolerance  =  1  x  10~8,  §p  =  90° 
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Table  1:  Computation  Time  (3.0GHz  Pentium  IV,  512MB  RAM) 


Error  Tolerance 

HPBW 

CPU  Time 

1  x  10-4 

3.38° 

18  sec 

1  x  10-5 

3.38° 

21  sec 

1  x  10^6 

3.38° 

49sec 

1  x  10-7 

3.38° 

2minl  lsec 

1  x  10-8 

3.38° 

5min22sec 

for  the  half  power  beam  width  for  the  different  error  tolerances.  Notice  that  the  HPBW  is 
unchanged  for  the  different  error  tolerances,  this  is  because  the  PO/Surface  Current  Method 
is  quite  accurate  for  the  main  beam  and  the  first  couple  of  side  lobes.  However,  if  we  wish 
to  use  this  technique  to  calculate  the  further  out  sidelobes  we  require  a  much  smaller  error 
tolerance,  which  is  clearly  depicted  here.  In  Figures  10  and  11  we  see  the  E  and  H-plane 
radiation  patterns  respectively.  Notice  that  the  cross-polarization  field  does  not  appear  in 
these  figures,  recall  that  the  feed  model  we  are  using  is  a  balanced  feed  where  we  have  a 
very  low  cross-polarized  field  level.  Hence  the  reflector  radiation  pattern  will  have  a  very 
low  cross-pol  level. 


Figure  10:  Radiation  Pattern  for  Reflector  #1 ,  4>p  =  0°,  MATLAB  Numerical  Integration 
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Relative  Power  Level  (dB) 
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,80 1 - ' - 1 - 1 - 1 - 1 - 

-30  -20  -10  0  10  20  30 


Angle  off  Broadside  -  0^  (degrees) 

Figure  11:  Radiation  Pattern  for  Reflector  #1 ,  ()>,,  =  90°,  MATLAB  Numerical  Integration 
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3.2.2  Reflector  #2 


We  will  now  look  at  the  results  of  this  model  with  a  displace  feed.  The  feed  is  being 
displaced  along  the  x-axis  or  alternatively  in  the  E-plane  ((^  =  0°).  We  would  then  expect 
to  the  main  beam  of  the  reflector  to  have  a  squint  in  the  E-plane.  This  amount  of  squint  can 
be  quantified  with  the  beam  deviation  factor  [5]  (BDF): 


eB  =  1+0.36  [4g]  2 

0f  i+^r2 


(60) 


where  0B  is  the  reflector  main  beam  scan  angle  and  0/.-  is  the  feed  displacement  angle.  For 
our  present  reflector  configuration  we  would  have  a  BDF  =  0.943  which  would  yield  a 
main  beam  squint  of  0.844°.  Figures  12  and  13  show  the  E  and  H-plane  radiation  patterns 
for  this  reflector  configuration  respectively.  We  can  see  from  Figure  12  that  there  is  a  beam 


Figure  12:  Radiation  Pattern  for  Reflector  #2,  (j)p  =  0°,  MATLAB  Numerical  Integration 


squint,  however  on  this  scale  it  is  not  apparent  on  the  exact  amount  of  squint.  Figure  14 
is  a  zoomed-in  view  of  the  main  beam  peak,  we  can  see  here  that  the  amount  of  beam 
squint  predicted  using  equation  (60)  is  quite  accurate,  the  model  predicted  a  main  beam  at 
—0.84°  off  broadside.  We  can  also  see  in  Figure  13  that  the  cross-pol  level  has  increased. 
This  is  expected  because  the  of  the  feed  displacement.  We  can  however  see  that  there 
is  a  significant  difference  between  the  computed  cross-pol  level  and  what  is  predicted  in 
GRASP. 
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Figure  13:  Radiation  Pattern  for  Reflector  # 2 ,  <\>p  =  90°,  MATLAB  Numerical  Integration 


Figure  14:  Radiation  Pattern  Beam  Squint  (Reflector  #2),  MATLAB  Numerical  Integration 
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3.2.3  Reflector  #3 


We  will  now  look  at  the  final  reflector  configuration  where  we  have  a  feed  displacement 
along  the  y-axis.  As  opposed  to  the  previous  case,  we  will  now  see  a  main  beam  squint  in 
the  H-plane  pattern.  As  before  the  BDF  remains  the  same  ( BDF  =  0.943),  however  with 
the  feed  displaced  by  X/2  we  should  see  a  main  beam  squint  of  1.69°.  In  Figures  15  and  16 
show  the  E  and  Fl-plane  patterns  respectively.  Figure  17  is  a  zoomed  in  view  of  the  main 


-30  -20  -10  0  10  20  30 


Angle  off  Broadside  -  0^  (degrees) 

Figure  15:  Radiation  Pattern  for  Reflector  #3,  MATLAB  Numerical  Integration,  Error 
Tolerance  =  1  x  10~8,  (j>p  =  0° 


beam  peak,  we  can  see  here  that  the  amount  of  beam  squint  predicted  is  accurate,  the  model 
predicted  a  main  beam  at  1.69°  off  broadside.  It  should  also  be  noted  that  because  we  have 
a  feed  displacement  the  integrand  of  the  radiation  integral  contains  additional  amplitude 
and  phase  terms.  Therefore  the  dblquad  routine  in  MATLAB  does  not  converge  as  quickly 
as  it  does  when  there  is  no  feed  displacement.  As  a  result  of  this  the  patterns  shown  in 
Figures  15  and  16  took  approximately  16  minutes  of  CPU  time  to  compute.  We  should  also 
notice  that  the  cross-pol  levels  are  ’’noisy”  in  Figures  15  and  16.  To  achieve  more  accurate 
cross-pol  patterns  we  would  have  to  decrease  the  error  tolerance.  Figures  18  and  19  show 
the  principal  plane  patterns  computed  using  an  error  tolerance  of  1  x  10  9.  We  can  see  that 
the  co  and  cross-pol  patterns  arc  predicted  accurately.  However  because  of  the  decreased 
error  tolerance  the  CPU  time  required  to  compute  each  of  these  patterns  was  well  over  one 
hour. 
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Figure  16:  Radiation  Pattern  for  Reflector  #3,  MATLAB  Numerical  Integration,  Error 
Tolerance  =  1  x  1CT8,  §p  =  90° 


Figure  17:  Radiation  Pattern  Beam  Squint  (Reflector  #3),  MATLAB  Numerical  Integration 
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Figure  18:  Radiation  Pattern  for  Reflector  #3,  MATLAB  Numerical  Integration,  Error 
Tolerance  =  1  x  10~9,  §p  =  0° 


Computed  Co-pol 
Grasp  Co-pol 


Angle  off  Broadside  -  0^  (degrees) 


Figure  19:  Radiation  Pattern  for  Reflector  #3,  MATLAB  Numerical  Integration,  Error 
Tolerance  =  1  x  10~9,  §p  =  90° 
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3.3  Ludwig’s  Method 


As  stated  previously  this  method  divides  the  integration  surface  into  a  grid  (M  x  A ).  Each 
cell  (incremental  surface  area)  of  this  grid  has  a  dimension  on  the  order  of  X  x  X.  The 
amplitude  and  phase  of  the  integrand  of  the  radiation  integral  is  then  approximated  by  a 
plane  over  this  incremental  surface  area.  If  we  select  the  size  of  this  incremental  area  such 
that  its  dimensions  are  on  the  order  of  one  wavelength,  we  should  have  accurate  results.  To 
do  this  we  need  to  determine  A0”!  and  A<|)”.  Consider  Figure  20.  Figure  20  (a)  shows  the 


Figure  20:  Reflector  Geometry  (a)  Projected  x-z  plane  (b)  Projected  x-y  plane 


projected  x-z  plane.  If  1  is  small  enough  we  can  safely  approximate  that  r\  =  m.  We  can 
then  determine  A0f  by  setting  /  =  X  and  using  the  formula  for  an  arc  length,  given  as 

/  =  rA0"!  (61) 

and 

2  F 

r=TT^e:  <62) 

0O  was  determined  in  Section  3.1.  We  can  then  determine  that  A0"!  ~  3.4708°,  this  gives 
that  M  =  10.  Figure  20  (b)  shows  the  projected  x-y  plane.  A(|)"  can  be  determined  in  much 
the  same  way  as  A0”'.  If  we  set  s  =  A,  and  using  the  equation  for  an  arc  length  and  solving 
for  A(j)" ,  we  get  A([>"  =  5.73°.  In  order  to  have  around  number  for  A  we  select  A  =  72  which 
gives  Acf)”  =  5°  and  an  arc  length  of  s  =  0.873A  Now  that  M  and  A  are  determined  we  can 
proceed  to  the  analysis  of  the  three  reflector  configurations.  As  stated  in  section  2.2  there  is 
a  need  for  large  amounts  of  memory  for  the  storage  of  coefficients  that  are  independent  of 
the  output  observation  angle.  With  this  in  mind  it  would  be  advantageous  for  this  routine  to 
be  coded  as  a  function  call  in  MATFAB  with  the  feed  location  and  the  pattern  cut  angle  (<j)p) 
as  input  variables.  Having  this  routine  coded  as  a  function  will  significantly  increase  the 
computational  efficiency  because  all  the  variables  of  the  script,  including  the  large  variables 
containing  the  coefficient  information,  are  local  to  the  function.  The  fact  that  these  large 
variables  are  not  part  of  the  MATFAB  workspace  will  decrease  the  time  required  to  run  this 
script. 
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3.3.1  Reflector  #1 


Figures  21  and  22  show  the  E  and  H-plane  patterns  respectively.  Both  these  patterns  were 
computed  in  only  7  seconds.  Note  that  both  these  patterns  closely  match  the  patterns  com¬ 


puted  with  GRASP.  We  should  also  notice  in  Figure  21  that  the  computed  cross-pol  level 
can  be  seen  at  approximately  -80dB  even  though  we  have  a  balanced  feed  and  should  have 
zero  cross-pol.  The  appearance  of  the  cross-pol  pattern  here  is  most  likely  due  to  numerical 
errors  introduced  by  approximating  the  amplitude  and  phase  terms  of  the  integrand  of  the 
radiation  integral  by  a  first  degree  polynomial  (i.e.  plane)  over  the  incremental  surface  area. 
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Relative  Power  Level  (dB) 
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Figure  22:  Radiation  Pattern  for  Reflector  #1,  §p  =  90°,  Ludwig  Method 
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3.3.2  Reflector  #2 


We  will  now  look  at  how  well  the  Ludwig  method  works  with  a  displaced  feed.  Figures  23 
and  24  show  the  E  and  Fl-plane  patterns  respectively  with  a  displaced  feed.  Once  again 
even  though  we  are  introducing  a  more  complicated  integrand,  both  these  patterns  were 
computed  in  7  seconds.  Again  with  the  displaced  feed  we  should  see  a  squint  in  the  main 


Computed  Co-pol 
Computed  X-pol 
Grasp  Co-pol 


Angle  off  Broadside  -  0p  (degrees) 


Figure  23:  Radiation  Pattern  for  Reflector  #2,  (f>„  =  0°,  Ludwig  Method 


beam  of  0.84°.  Figure  25  shows  a  zoomed  in  view  of  the  squinted  main  beam.  We  can  see 
that  the  Ludwig  method  accurately  models  the  squint  in  the  main  beam. 
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Figure  24:  Radiation  Pattern  for  Reflector  #2,  <\>p  =  90°,  Ludwig  Method 
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3.3.3  Reflector  #3 


We  will  now  look  at  the  final  reflector  configuration.  Figures  26  and  27  show  the  E  and 
Fl-plane  patterns  respectively.  As  it  was  the  case  for  the  other  two  configurations  each  of 
these  pattern  cuts  were  produced  in  only  7  seconds.  Notice  here  that  cross-pol  patterns  are 


modelled  fairly  accurately.  The  numerical  errors  introduced  because  of  the  approximation 
made  by  the  Ludwig  method  does  not  seem  to  be  a  significant  issue  when  an  actual  cross- 
pol  pattern  is  modelled.  Again  we  see  a  main  beam  squint.  Figure  28  shows  an  zoomed 
in  view  of  the  main  beam  and  we  see  again  that  the  Ludwig  method  accurately  models  the 
squint  in  the  main  beam. 
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Figure  28:  Radiation  Pattern  Beam  Squint  (Reflector  #3),  Ludwig  Method 
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4  Concluding  Remarks 


We  can  see  from  the  results  presented  in  the  previous  section  that  both  of  the  methods 
discussed  (PO/Surface  Current  Method  and  Ludwig’s  Method)  give  accurate  results  when 
compared  to  patterns  computed  using  commercially  available  reflector  modelling  software. 
Where  these  two  methods  differ  is  in  the  amount  of  time  required  to  computed  these  re¬ 
sults.  We  have  seen  that  with  the  PO/SC  Method,  the  integration  of  the  surface  currents 
induced  on  the  reflector  surface  due  to  the  feed  is  performed  using  a  numerical  integration 
routine  built  into  MATLAB.  We  have  shown  that,  with  a  laterally  displaced  feed,  where  the 
integrand  of  the  radiation  integral  increases  in  complexity,  the  evaluation  of  one  principal 
plane  pattern  can  take  up  to  one  hour  to  compute.  Incorporating  the  PO/SC  method  into  a 
seeker  model  would  therefore  be  impractical. 

As  an  alternative  to  the  PO/SC  Method  we  looked  at  Ludwig’s  Method  for  the  analysis 
of  reflector  patterns.  We  have  shown  that  this  method,  regardless  of  the  lateral  displace¬ 
ment  of  the  feed,  can  accurately  compute  the  principal  plane  pattern  in  under  10  seconds. 
It  was  also  shown  that  this  method  did  introduce  some  numerical  errors,  specifically  in  the 
cross  polarization  pattern,  due  to  the  approximations  made.  Further  investigation  will  be 
required  to  see  if  these  numerical  errors  cause  significant  inaccuracies  and  if  so  how  to 
eliminate  or  reduce  these  errors  so  that  they  do  not  cause  major  issues. 

Future  work  will  involve  applying  these  modelling  techniques  to  reflector  antennas  which 
have  a  compound  feed  for  use  in  monopluse  systems.  Finally  a  monopulse  antenna  model, 
based  on  the  Ludwig  method,  could  be  incorporated  into  a  missile  seeker  model. 
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Annex  A 
MATLAB  Code 


clear  clc 

%**  *******************************************  ***************  ********** 
0 

Q, 

0 

%File  name:  Ref lector_Rad_Pattern_PO .m 

0, 

0 

%This  MATLAB  routine  will  compute  the  radiation  pattern  of  a  parabolic 

%ref lector  using  the  PO/Surface  Current  Method. 

0, 

0 

o 

%Global  variables 

global  EO  UO  ZO  C  FREQ  LAM  W  KO  D  F 

%Start  time  to  determine  run  time 
S_time  =  clock; 

%Defines  constants 

EO  =  8 . 8542e-12;  UO  =  4*pi*le-7;  ZO  =  sqrt (UO/EO) ;  C  = 

(sqrt (EO*UO) ) "-1; 

%Model  electrical  parameters 

FREQ  =  30e9;  LAM  =  C/FREQ;  W  =  2*pi*FREQ;  KO  =  W*sqrt (UO*EO) ; 

%Reflector  physical  parameters 
D  =  20*LAM;  FoD  =  0.8;  F  =  FoD*D; 

%Far-Field  distance 
R  =  (  (2*D~2) /LAM) ; 

%Determines  the  subtended  angle  of  the  parabaloid 
theta_o_s  =  2*atan2  (D, 4*F ) ; 

%Vector  which  defines  the  feed  position 
F_v  =  [ 0 ; 0 ; 0 ] ; 

%Feed  raised  cosine  factor 
qe  =  6.5;  qh  =  6.5; 

%Radiation  pattern  cut 

phi_p_deg  =  0;  phi_p  =  rad (phi_p_deg) ; 
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o\°  o\° 


%Determines  repeated  cos  and  sin  functions  of  phi_p 

cpp  =  cos (phi _ p ) ;  spp  =  sin(phi_p);  cpp2  =  cpp~2;  spp2  =  spp"2; 

%Constant  definition 

const  =  (-j*KO*ZO*exp(-j*KO*R) )/ (4*pi*R) ; 

%Index  initiallization 
ii  =  1; 

For  loop  which  will  calculate  the  radiation  pattern  as  a  function  of  the 
angle  off  broadside  (theta_p) 
for  theta_p_deg  =  -30:0.2:30 
theta_p  =  rad (theta_p_deg) ; 

%Determines  repeated  cos  and  sin  functions  of  theta_p 

stp  =  sin (theta_p) ; 

ctp  =  cos (theta_p) ; 

stp2  =  stp~2; 

ctp2  =  ctp~2; 

%Function  calls  (dblquad)  that  will  numerically  evaluate  the  integrand 
%of  the  radiation  integral  for  the  3  cartesian  coordinates. 

Fx  =  dblquad (@Rad_int_x_rev4, 0, 2*pi, 0, theta_o_s, le-8, @quad, F_v, ... 
qe,  qh,  stp, ctp, spp, cpp) ; 

Fy  =  dblquad (@Rad_int_y_rev4, 0, 2*pi, 0, theta_o_s, le-8, @quad, F_v, ... 
qe,  qh,  stp, ctp, spp, cpp) ; 

Fz  =  dblquad (@Rad_int_z_rev4, 0, 2*pi, 0, theta_o_s, le-8, @quad, F_v, ... 
qe,  qh,  stp, ctp, spp, cpp) ; 

%Computes  the  cartesian  components  of  the  far-zone  electic  field 
% (radiation  pattern) 

Ex  =  const* (Fx  -  Fx*stp2*cpp2  +  Fy*stp2*spp*cpp  +  Fz*ctp*stp*cpp) ; 

Ey  =  const* (-Fy  -  Fx*stp2*cpp*spp  +  Fy*stp2*spp2  +  Fz*ctp*stp*spp) ; 

Ez  =  const*  (-Fz  -  Fx*stp*cpp*ctp  +  Fy*stp*spp*ctp  +  Fz*ctp2) ; 

%Computes  the  co  and  x-pol  patterns  bassed  on  the  cartessian  components 
%determined  previously 

E_co_p(ii)  =  -(1  -  ctp) *spp*cpp*Ex  +  (1  -  spp2*  (1  -  ctp))*Ey... 

-  stp*spp*Ez; 

E_x_p(ii)  =  +  (1  -  cpp2*(l  -  ctp) ) *Ex  -  (1  -  ctp) *spp*cpp*Ey . . . 

-  stp*cpp*Ez; 

%Increment  index 
ii  =  ii  +  1; 

end 


DRDC  Ottawa  TM  2005-098 


33 


%Defines  theta_p  to  plot  pattern 
theta  =  -30:0.2:30; 

%Determines  the  co  and  x-pol  normalized  patterns  in  dB 
m_E_co  =  abs(E_co_p);  m_E_x  =  abs (E_x_p) ;  peaks  =  [max(m_E_co) 
max(m_E_x)];  max_val  =  max (peaks); 

n_E_co  =  m_E_co/max_val;  n_E_x  =  m_E_x/max_val; 

db_E_co  =  20*logl0 (n_E_co) ;  db_E_x  =  20*logl0 (n_E_x) ; 

%Plots  normalized  radiation  patterns 

figure  plot (theta, db_E_co, ' -k' LineWidth' , 2)  hold  on 

plot (theta, db_E_x, ' -r' , ' LineWidth' , 2) 

%Defines  figure  parameters 

ylim  ([-80  0])  xlabel ('Angle  off  Broadside  -  \theta_p 
(degrees) ' , ' FontName' ,  .  .  . 

'Times  New  Roman' ,' FontSize' , 14 ,' FontWeight' ,' bold' ) 
ylabel (' Relative  Power  Level  (dB) ',' FontName' . 

'Times  New  Roman' ,' FontSize' , 14 ,' FontWeight' ,' bold' ) 
grid  on 

%Routine  that  stops  timer  and  converts  elapsed  time  to  h:m:s  format 

t ime ( S_t ime ) 
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function  [E_temp]  = 

Rad_int_x_rev4 (phi_s, theta_s, F_v, qe, qh, stp, ctp, spp, cpp) 

o 

o. 

o 

%File  name:  Rad_int_x_rev4 .m 

0. 

0 

%This  function  will  define  the  x-component  of  the  integrand  which  will  be 
%numerically  integrated  using  dblquad  in  Ref lector_Rad_Pattern_PO .m 

Q. 

0 

%**  ****************************************************************  ****** 
0 

%Global  variables 

global  EO  UO  ZO  C  FREQ  LAM  W  KO  D  F 
%Defines  repeated  sinusoidal  terms 

sts  =  sin (theta_s) ;  cts  =  cos (theta_s) ;  sps  =  sin(phi_s);  cps  = 
cos(phi_s);  sts2  =  sin (theta_s/2) ;  cts2  =  cos (theta_s/2) ;  tts2  = 
tan (theta_s/2) ; 

%Defines  rho  vector  from  origin  to  reflector  surface 

r_s  =  (2*F)./(1  +  cts);  r_sx  =  r_s . *sts*cps;  r_sy  =  r_s . *sts*sps; 

r_sz  =  r_s.*cts; 

%Defines  angles  of  primed  coor  sys . 

terml  =  sqrt ( (r_s . *sts*cps  -  F_v(l))."2  +  (r_s . *sts*sps  - 
F_v (2) ) . "2) ;  term2  =  r_s.*cts;  theta_s_p  =  atan2 (terml, term2) ; 

term3  =  (r_s . *sts*sps  -  F_v(2));  term4  =  (r_s . *sts*cps  -  F_v(l)); 
phi_s_p  =  atan2 (term3, term4) ; 

%Defines  repeated  sinusoidal  terms 

stsp  =  sin (theta_s_p) ;  ctsp  =  cos (theta_s_p) ;  spsp  =  sin (phi_s_p) ; 
cpsp  =  cos (phi_s_p) ; 

%Defines  the  magnitude  of  rho  prime  vector  which  extends  from  the  feed  to 

%the  surface  of  the  reflector 

r_s_p  =  sqrt (r_s . "2  +  F_v(l)"2  +  F_v(2)~2  ... 

-  2*r_s . *sts . * (F_v (1) *cps  +  F_v (2) *sps) ) ; 
const  =  (2*exp (- j*KO*r_s_p) ) / (r_s_p*ZO) ; 

%Defines  the  feed  pattern  coefficients 
Ce  =  (ctsp).'qe;  Ch  =  (ctsp).'qh; 

%Defines  the  x  component  of  the  surface  current  induced  on  the  reflector 
%surface 
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Js_x  =  const* (-Ch . *cpsp . *stsp . *sts2 . *sps  ..  . 

-Ch . *cpsp . *ctsp . *spsp . *cts2 . . . 

+Ce . *spsp . *cpsp . *cts2 ) ; 

%Defines  exponential  term  of  the  integrand 

tempi  =  r_s . * (sts*cps*stp*cpp  -  sts*sps*stp*spp  -  cts*ctp) ;  expo  = 
exp ( j*KO*templ ) ; 

%Defines  differential  surface  term 
ds  =  r_s . "2 . *sts . *sec (theta_s/2) ; 

%Defines  the  integrand  to  be  numerically  integrated 
E_temp  =  Js_x . *expo . *ds; 
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function  [E_temp]  = 

Rad_int_y_rev4 (phi_s, theta_s, F_v, qe, qh, stp, ctp, spp, cpp) 

^************************************************************************ 

o 

o. 

o 

%File  name:  Rad_int_y_rev4 .m 

0. 

0 

%This  function  will  define  the  y-component  of  the  integrand  which  will  be 
%numerically  integrated  using  dblquad  in  Ref lector_Rad_Pattern_PO .m 

Q. 

0 

%**  ****************************************************************  ****** 
0 

%Global  variables 

global  EO  UO  ZO  C  FREQ  LAM  W  KO  D  F 
%Defines  repeated  sinusoidal  terms 

sts  =  sin (theta_s) ;  cts  =  cos (theta_s) ;  sps  =  sin(phi_s);  cps  = 
cos(phi_s);  sts2  =  sin (theta_s/2) ;  cts2  =  cos (theta_s/2) ;  tts2  = 
tan (theta_s/2) ; 

%Defines  rho  vector  from  origin  to  reflector  surface 

r_s  =  (2*F)./(1  +  cts);  r_sx  =  r_s . *sts*cps;  r_sy  =  r_s . *sts*sps; 

r_sz  =  r_s.*cts; 

%Defines  angles  of  primed  coor  sys . 

terml  =  sqrt ( (r_s . *sts*cps  -  F_v(l))."2  +  (r_s . *sts*sps  - 
F_v (2) ) . "2) ;  term2  =  r_s.*cts;  theta_s_p  =  atan2 (terml, term2) ; 

term3  =  (r_s . *sts*sps  -  F_v(2));  term4  =  (r_s . *sts*cps  -  F_v(l)); 
phi_s_p  =  atan2 (term3, term4) ; 

%Defines  repeated  sinusoidal  terms 

stsp  =  sin (theta_s_p) ;  ctsp  =  cos (theta_s_p) ;  spsp  =  sin (phi_s_p) ; 
cpsp  =  cos (phi_s_p) ; 

%Defines  the  magnitude  of  rho  prime  vector  which  extends  from  the  feed  to 

%the  surface  of  the  reflector 

r_s_p  =  sqrt (r_s . "2  +  F_v(l)"2  +  F_v(2)~2  ... 

-  2*r_s . *sts . * (F_v (1) *cps  +  F_v (2) *sps) ) ; 
const  =  (2*exp (- j*KO*r_s_p) ) / (r_s_p*ZO) ; 

%Defines  the  feed  pattern  coefficients 
Ce  =  (ctsp).'qe;  Ch  =  (ctsp).'qh; 

%Defines  the  y  component  of  the  surface  current  induced  on  the  reflector 
%surface 
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Js_y  =  const* (Ch . * (cpsp) . "2 . *ctsp . *cts2 . . . 

+  Ce . * (spsp) . "2 . *cts2  .  .  . 

+  Ch . *cpsp .  *stsp .  *sts2  .  *cps) ; 

%Defines  exponential  term  of  the  integrand 

tempi  =  r_s . * (sts*cps*stp*cpp  -  sts*sps*stp*spp  -  cts*ctp) ;  expo  = 
exp ( j*KO*templ ) ; 

%Defines  differential  surface  term 
ds  =  r_s . *2 . *sts . *sec (theta_s/2) ; 

%Defines  the  integrand  to  be  numerically  integrated 
E_temp  =  Js_y . *expo . *ds; 
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function  [E_temp]  = 

Rad_int_z_rev4 (phi_s, theta_s, F_v, qe, qh, stp, ctp, spp, cpp) 

o 

o. 

o 

%File  name:  Rad_int_z_rev4 .m 

0. 

0 

%This  function  will  define  the  z-component  of  the  integrand  which  will  be 
%numerically  integrated  using  dblquad  in  Ref lector_Rad_Pattern_PO .m 

Q. 

0 

%**  ****************************************************************  ****** 
0 

%Global  variables 

global  EO  UO  ZO  C  FREQ  LAM  W  KO  D  F 
%Defines  repeated  sinusoidal  terms 

sts  =  sin (theta_s) ;  cts  =  cos (theta_s) ;  sps  =  sin(phi_s);  cps  = 
cos(phi_s);  sts2  =  sin (theta_s/2) ;  cts2  =  cos (theta_s/2) ;  tts2  = 
tan (theta_s/2) ; 

%Defines  rho  vector  from  origin  to  reflector  surface 

r_s  =  (2*F)./(1  +  cts);  r_sx  =  r_s . *sts*cps;  r_sy  =  r_s . *sts*sps; 

r_sz  =  r_s.*cts; 

%Defines  angles  of  primed  coor  sys . 

terml  =  sqrt ( (r_s . *sts*cps  -  F_v(l))."2  +  (r_s . *sts*sps  - 
F_v (2) ) . "2) ;  term2  =  r_s.*cts;  theta_s_p  =  atan2 (terml, term2) ; 

term3  =  (r_s . *sts*sps  -  F_v(2));  term4  =  (r_s . *sts*cps  -  F_v(l)); 
phi_s_p  =  atan2 (term3, term4) ; 

%Defines  repeated  sinusoidal  terms 

stsp  =  sin (theta_s_p) ;  ctsp  =  cos (theta_s_p) ;  spsp  =  sin (phi_s_p) ; 
cpsp  =  cos (phi_s_p) ; 

%Defines  the  magnitude  of  rho  prime  vector  which  extends  from  the  feed  to 

%the  surface  of  the  reflector 

r_s_p  =  sqrt (r_s . "2  +  F_v(l)"2  +  F_v(2)~2  ... 

-  2*r_s . *sts . * (F_v (1) *cps  +  F_v (2) *sps) ) ; 
const  =  (2*exp (- j*KO*r_s_p) ) / (r_s_p*ZO) ; 

%Defines  the  feed  pattern  coefficients 
Ce  =  (ctsp).'qe;  Ch  =  (ctsp).'qh; 

%Defines  the  z  component  of  the  surface  current  induced  on  the  reflector 
%surface 
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Js_z  =  const* (Ch . *cpsp . *ctsp . *spsp . *sts2 . *cps .. . 

-Ce . *spsp . *cpsp . *sts2 . *cps  .  .  . 

-Ch . * (cpsp) . "2 . *ctsp . *sts2 . *sps  .  .  . 

-Ce . * (spsp) . "2 . *sts2 . *sps) ; 

%Defines  exponential  term  of  the  integrand 

tempi  =  r_s . * (sts*cps*stp*cpp  -  sts*sps*stp*spp  -  cts*ctp) ;  expo  = 
exp ( j*KO*templ ) ; 

%Defines  differential  surface  term 
ds  =  r_s . "2 . *sts . *sec (theta_s/2) ; 

%Defines  the  integrand  to  be  numerically  integrated 
E_temp  =  Js_z . *expo . *ds; 
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function  [E_co_p, E_x_p]  = 

Ref lector_Rad_Pattern_Ludwig_FC (F_v, phi_p_deg) 

o 

o. 

o 

%File  name:  Ref lector_Rad_Pattern_Ludwig_FC .m 

0. 

0 

%This  MATLAB  function  will  compute  the  radiation  pattern  of  a  parabolic 

%reflector  using  the  Ludwig  Method. 

0, 

0 

%-** ****************************************************************  ***** 
global  EO  UO  ZO  C  FREQ  LAM  W  KO  D  F 

%Far-Field 
R  =  (2*D"2) /LAM; 

%Determines  the  subtended  angle  of  the  parabaloid 
theta_o_s  =  2*atan2 (D, 4*F ) ; 

%Feed  pattern  raised  cosine  coefficients 
qe  =  6.5;  qh  =  6.5; 

%Defines  the  grid  on  the  reflecotr  surface  which  defines  the  incremental 
%surface  area 

theta_seg  =  11;  del_theta  =  theta_o_s/ (theta_seg-l) ;  theta_t  =  0; 
for  mm  =  1 : 1 :theta_seg 

theta_s_i (mm)  =  theta_t; 

theta_t  =  theta_s_i (mm)  +  del_theta; 

end 

phi_seg  =  73;  del_phi  =  2*pi/ (phi_seg-l) ;  phi_t  =  0;  for  nn  = 

1 : 1 :phi_seg 

phi_s_i(nn)  =  phi_t; 

phi_t  =  phi_s_i(nn)  +  del_phi; 

end 

%Computes  repeated  sinusoidal  terms 

stsi  =  sin (theta_s_i) ;  ctsi  =  cos (theta_s_i) ;  spsi  =  sin (phi_s_i) ; 
cpsi  =  cos (phi_s_i) ; 

%Computes  length  from  the  origine  to  reflector  surface 
r_s  =  (2*F)./(1  +  cos (theta_s_i) ) ; 

%Radiation  pattern  cut  in  radians 
phi_p  =  rad (phi_p_deg) ; 
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%Computes  repeated  sinusoidal  terms 

spp  =  sin(phi_p);  cpp  =  cos  (phi_p) ;  cpp2  =  cpp~2;  spp2  =  spp"2; 

%Preallocation  of  storage  vectors 
r_s_pl  =  zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fxmn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fymn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fzmn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mnl  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mn2  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mn3  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ; 

r_s_p2  =  zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fxmpn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fympn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fzmpn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mpnl  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mpn2  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mpn3  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ; 

r_s_p3  =  zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fxmnp  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fymnp  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fzmnp  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mnpl  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mnp2  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mnp3  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ; 

r_s_p4  =  zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fxmpnp  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fympnp  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  Fzmpnp  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mpnpl  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mpnp2  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  temp_mpnp3  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ; 

ax_mn  =  zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  bx_mn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  cx_mn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  ay_mn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  by_mn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  cy_mn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  az_mn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  bz_mn  = 
zeros ( (theta_seg-l) * (phi_seg-l) , 1) ;  cz_mn  = 
zeros ( (theta_seg-l) * (phi_seg-l) ,  1) ; 
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%Index  initiallization 

j  j  =  i  ; 


%This  routine  will  calculate  all  the  terms  that  are  independant  on  the 
%output  pattern  angle  and  stores  them  in  vectors.  Namely  it  will  compute 
%the  amplitude  term  at  all  the  corners  of  the  incremental  surface  area  for 
%the  three  cartesian  components  it  also  computes  the  coefficients  of  the 
%linear  function  which  approximates  the  amplitude  over  the  incremental 
%surface  area, 
for  mm  =  1 : 1 :theta_seg-l 
for  nn  =  1 : 1 :phi_seg-l 

[ r_s_pl ( j  j ) , Fxmn ( j  j ) , Fymn ( j  j ) , Fzmn ( j  j ) , temp_mnl ( j  j ) , . . . 

temp_mn2 ( j j ) , temp_mn3 ( j j ) ]  =  Ludwig_terms (mm, nn, theta_s_i, . . . 
st si,  ct si,  spsi,  cpsi,  F_v,  r_s,  qe,  qh)  ; 

[ r_s_p2 ( j  j ) , Fxmpn ( j  j ) , Fympn ( j  j ) , Fzmpn ( j  j ) , temp_mpnl ( j  j ) , . . . 
temp_mpn2 ( j j ) , temp_mpn3 ( j j ) ]  =  Ludwig_terms ( (mm+1) , nn, . . . 
theta_s_i, stsi, ctsi, spsi, cpsi, F_v, r_s, qe, qh) ; 

[ r_s_p3 ( j  j ) , Fxmnp ( j  j ) , Fymnp ( j  j ) , Fzmnp ( j  j ) , temp_mnpl ( j  j ) , . . . 
temp_mnp2 ( j j) ,temp_mnp3 ( j j) ]  =  Ludwig_terms (mm, (nn+1) , . . . 
theta_s_i, stsi, ctsi, spsi, cpsi, F_v, r_s, qe, qh) ; 

[ r_s_p 4 ( j  j ) , Fxmpnp ( j  j ) , Fympnp ( j  j ) , Fzmpnp ( j  j ) , temp_mpnpl ( j  j ) , . . . 
temp_mpnp2 ( j j) ,temp_mpnp3 ( j j) ]  =  Ludwig_terms ( (mm+1) , . . . 

(nn+1) , theta_s_i, stsi, ctsi, spsi, cpsi, F_v, r_s, qe, qh) ; 

ax_mn(jj)  =  (1/4) * (3*Fxmn ( j j)  -  Fxmpnp ( j j ).. . 

+  Fxmpn (jj)  +  Fxmnp(jj)); 
bx_mn(jj)  =  (1/ (2*del_theta) )* (Fxmpn ( j j ).. . 

-  Fxmn(jj)  +  Fxmpnp (jj)  -  Fxmnp(jj)); 
cx_mn(jj)  =  (1/ (2*del_phi) )* (Fxmnp (jj) .. . 

-  Fxmn(jj)  +  Fxmpnp (jj)  -  Fxmpn(jj)); 

ay_mn(jj)  =  (1/4) * (3*Fymn ( j j)  -  Fympnp ( j j ).. . 

+  Fympn (jj)  +  Fymnp (jj) ) ; 
by_mn(jj)  =  (1/ (2*del_theta) )* (Fympn ( j j ).. . 

-  Fymn ( j  j )  +  Fympnp (jj)  -  Fymnp(jj)); 
cy_mn(jj)  =  (1/ (2*del_phi) )* (Fymnp (jj) .. . 

-  Fymn ( j  j )  +  Fympnp (jj)  -  Fympn(jj)); 

az_mn(jj)  =  (1/4) * (3*Fzmn ( j j)  -  Fzmpnp ( j j ).. . 

+  Fzmpn (jj)  +  Fzmnp(jj)); 
bz_mn(jj)  =  (1/ (2*del_theta) )* (Fzmpn ( j j ).. . 

-  Fzmn(jj)  +  Fzmpnp (jj)  -  Fzmnp(jj)); 
cz_mn(jj)  =  (1/ (2*del_phi) )* (Fzmnp ( j j ).. . 

-  Fzmn(jj)  +  Fzmpnp (jj)  -  Fzmpn(jj)); 
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jj  =  jj  +i; 


end 

end 

%Defines  the  angles  off  broadside  for  which  the  radiation  pattern  is 
%computed 

theta  =  -30:0.2:30; 

%Preallocation  of  storage  vectors 

L  =  length (theta) ;  E_theta_p  =  zeros  ( L , 1 ) ;  E_phi_p  =  zeros (L,l); 

E_co_p  =  zeros  ( L, 1 ) ;  E_x_p  =  zeros (L,l); 

%Constant  definition 

const  =  (-j*KO*ZO*exp(-j*KO*R) )/ (4*pi*R) ; 

%Index  initiallization 
ii  =  1; 

%For  loop  which  will  calculate  the  radiation  pattern  as  a  function  of  the 
%angle  off  broadside  (theta_p) 
for  theta_p_deg  =  -30:0.2:30 
theta_p  =  rad (theta_p_deg) ; 

%Computes  repeated  sinusoidal  terms 

stp  =  sin (theta_p) ; 

ctp  =  cos (theta_p) ; 

stp2  =  stp~2; 

ctp2  =  ctp~2; 

%Variable  initialisation 
Fx  =  0; 

Fy  =  0; 

Fz  =  0; 

%Index  initiallization 

jj  =  i; 

%This  routine  computes  the  phase  term  of  the  integrand  on  the 
%integration  grid  which  defines  the  incremental  surface  area  it  also 
%computes  the  coefficients  linear  function  approximating  the  phase 
%onver  the  incremental  surface  area.  Finally  this  rouitne  also 
%computes  the  intgral  (closed  form)  over  the  incremental  surface  area 
%and  sums  all  these  contributions  for  the  3  cartesian  components, 
for  mm  =  1 : 1 :theta_seg-l 
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for  nn  =  1 : 1 :phi_seg-l 


gmn  =  temp_mnl ( j j ) *stp*cpp . . . 

-  temp_mn2 ( j j ) *stp*spp . . . 

-  temp_mn3 ( j j ) *ctp  -  r _ s _ p 1 ( j  j ) ; 

gmpn  =  temp_mpnl ( j j ) *stp*cpp .  .  . 

-  temp_mpn2 ( j j ) *stp*spp . . . 

-  temp_mpn3 ( j j ) *ctp  -  r_s_p2(jj); 

gmnp  =  temp_mnpl ( j j ) *stp*cpp .  .  . 

-  temp_mnp2 ( j j ) *stp*spp . .  . 

-  temp_mnp3 ( j j ) *ctp  -  r _ s _ p3 ( j  j ) ; 

gmpnp  =  temp_mpnpl ( j j ) *stp*cpp . . . 

-  temp_mpnp2 ( j j ) *stp*spp .  .  . 

-  temp_mpnp3 ( j j ) *ctp  -  r _ s _ p4 ( j  j ) ; 

alpha_mn  =  (l/4)*(3*gmn  -  gmpnp  +  gmpn  +  gmnp); 
beta_mn  =  (1/ (2*del_theta) ) * (gmpn  -  gmn  +  gmpnp  -  gmnp); 
xi_mn  =  (1/ (2*del_phi) ) * (gmnp  -  gmn  +  gmpnp  -  gmpn); 

terml  =  exp ( j*KO*alpha_mn) ; 
if  beta_mn  ==  0  |  xi_mn  ==  0 
term2  =  0; 
term3  =  0; 
term4  =  0; 

else 

term2  =  ( (exp ( j*KO*beta_mn*del_theta) -1) / ( j*KO*beta_mn) ) . 

* ( (exp ( j*KO*xi_mn*del_theta) -1) / ( j*K0*xi_mn) ) ; 
term3  =  ( ( (del_theta) / ( j*KO*beta_mn) ) . . . 

*exp ( j*KO*beta_mn*del_theta) . . . 

-  ( (exp ( j*KO*beta_mn*del_theta) -1) / ( j*KO*beta_mn) "2) ) 
* ( (exp ( j*KO*xi_mn*del_theta) -1) / ( j*K0*xi_mn) ) ; 

term4  =  ( (exp ( j*KO*beta_mn*del_theta) -1) / ( j*KO*beta_mn) ) . 
* ( ( (del_phi) / ( j*K0*xi_mn) ) *exp ( j*KO*xi_mn*del_phi) . . . 

-  ( (exp ( j*KO*xi_mn*del_phi) -1) / ( j*K0*xi_mn) "2) ) ; 

end 

del_Fx  =  terml* (ax_mn ( j j ) *term2  +  bx_mn ( j j ) *term3 . . . 

+  cx_mn ( j j ) *term4 ) ; 

del_Fy  =  terml* (ay_mn ( j j ) *term2  +  by_mn ( j j ) *term3 . . . 

+  cy_mn ( j j ) *term4 ) ; 

del_Fz  =  terml* (az_mn ( j j ) *term2  +  bz_mn ( j j ) *term3 . . . 

+  cz_mn ( j j ) *term4 ) ; 
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Fx  =  Fx  +  del_Fx; 
Fy  =  Fy  +  del_Fy; 
Fz  =  Fz  +  del_Fz; 

jj  =  jj  +  i; 


end 

end 

%Computes  the  cartesian  components  of  the  far-zone  electic  field 
% (radiation  pattern) 

Ex  =  const*  (Fx  -  Fx*stp2*cpp2  +  Fy*stp2*spp*cpp  +  Fz*ctp*stp*cpp) ; 

Ey  =  const*  (-Fy  -  Fx*stp2*cpp*spp  +  Fy*stp2*spp2  +  Fz*ctp*stp*spp) ; 

Ez  =  const*  (-Fz  -  Fx*stp*cpp*ctp  +  Fy*stp*spp*ctp  +  Fz*ctp2) ; 

%Computes  the  co  and  x-pol  patterns  bassed  on  the  cartessian  components 
%determined  previously 

E_co_p(ii)  =  -(1  -  ctp) *spp*cpp*Ex  +  (1  -  spp2*  (1  -  ctp))*Ey... 

-  stp*spp*Ez; 

E_x_p(ii)  =  +  (1  -  cpp2*(l  -  ctp) ) *Ex  -  (1  -  ctp) *spp*cpp*Ey . . . 

-  stp*cpp*Ez; 

%Increments  index 
ii  =  ii  +  1; 

end 


46 


DRDC  Ottawa  TM  2005-098 


function  [r_s_p, Fx, Fy, Fz, temp_l, temp_2, temp_3] . . . 

=  Ludwig_terms (mm, nn, theta_s_i, stsi, ctsi, spsi, cpsi, F_v, r_s, qe, qh) 

o 

o. 

o 

%File  name:  Ludwig_terms .m 

0. 

0 

%This  MATLAB  routine  will  compute  several  terms  used  in  the  Ludwig  Method 
%specif ically  the  terms  that  independant  on  the  output  pattern  angle 

%direction.  This  is  computed  for  differenc  m  and  n  indices. 

0, 

0 

o 

%Global  variables 

global  EO  UO  ZO  C  FREQ  LAM  W  KO  D  F 

%Defines  angles  of  primed  coor  sys .  (Feed  coor  sys . ) 
terml  =  sqrt ( (r_s (mm) *stsi (mm) *cpsi (nn)  -  F_v(l))~2  ... 

+  (r_s (mm) *stsi (mm) *spsi (nn)  -  F_v(2))"2); 
term2  =  r_s (mm) *ctsi (mm) ;  theta_s_p_i  =  atan2 (terml, term2) ; 

term3  =  (r_s (mm) *stsi (mm) *spsi (nn)  -  F_v(2));  term4  = 

(r_s (mm) *stsi (mm) *cpsi (nn)  -  F_v(l));  phi_s_p_i  = 
atan2 (term3, term4 ) ; 

%Defines  repeated  sinusoidal  terms 

stspi  =  sin (theta_s_p_i) ;  ctspi  =  cos (theta_s_p_i) ;  spspi  = 
sin (phi_s_p_i) ;  cpspi  =  cos (phi_s_p_i) ; 

%Defines  the  magnitude  of  rho  prime  vector  which  extends  from  the  feed  to 
%the  surface  of  the  reflector 

r_s_p  =  sqrt  (r_s  (mm)  "2  +  F_v(l)/'2  +  F_v(2)"2  ... 

-  2*r_s (mm) *stsi (mm) * (F_v (1) *cpsi (nn)  +  F_v (2) *spsi (nn) ) ) ; 

%Computes  the  amplitude  term  of  the  integrand  for  the  3  cartesian  coord. 
Fx  =  r_s (mm) "2*sec (theta_s_i (mm) /2) *stsi (mm) . . . 

* (-ctspi~qh*sin (theta_s_i (mm) /2) *spsi (nn) *stspi*cpspi . . . 

-ctspi"qh*cos (theta_s_i (mm) /2) *ctspi*cpspi*spspi . . . 

+ctspi'qe*cos (theta_s_i (mm) /2) *cpspi*spspi) * (2/ (ZO*r_s_p) ) ; 

Fy  =  r_s (mm) "2*sec (theta_s_i (mm) /2) *stsi (mm) . . . 

* (ctspi~qh*cpspi"2*ctspi*cos (theta_s_i (mm) /2) . . . 

+  ctspi~qe*spspi"2*cos (theta_s_i (mm) /2) ... 

+  ctspi"qh*cpspi*stspi*sin (theta_s_i (mm) /2) *cpsi (nn) ) * (2/ (ZO*r_s_p) ) ; 
Fz  =  r_s (mm) "2*sec (theta_s_i (mm) /2) *stsi (mm) . . . 
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* (ctspi~qh*cpspi*ctspi*spspi*sin (theta_s_i (mm) /2) *cpsi (nn) . . . 
-ctspi"qe*spspi*cpspi*sin (theta_s_i (mm) / 2 ) *cpsi (nn) . . . 

-ctspi'qh*cpspi "2 . *ctspi*sin (theta_s_i (mm) /2)*spsi(nn) . . . 
-ctspi'qe*spspi "2 . *sin (theta_s_i (mm) /2)*spsi(nn))*(2/ (ZO*r_s_p) ) ; 

%Computes  the  cartesian  components  of  the  point  on  the  reflector  where  the 

%amplitude  and  phase  terms  of  the  integrands  will  be  calculated  (i.e. 

%corners  of  the  incremental  surface  area. 

temp_l  =  r_s (mm) *stsi (mm) *cpsi (nn) ;  temp_2  = 

r_s (mm) *stsi (mm) *spsi (nn)  ;  temp_3  =  r_s (mm) *ctsi (mm) ; 
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